home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / r_check2du.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  4.8 KB  |  224 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3.  
  4. #include "fft.h"
  5. #include "constant.h"
  6.  
  7. /*                        */
  8. /*    Precision dependant declarations    */
  9. /*                        */
  10. #ifdef    DOUBLE
  11.  
  12. #define TOLERANCE   DTOLERANCE
  13. typedef        double    this_half;
  14. typedef        double    this_type;
  15.  
  16. #define        THIS_SUM    dsum_
  17. #define        THIS_GEN    dgen_
  18. #define        THIS_FTU    dft2du_
  19.  
  20. #define        THIS_FFTUI    dfft2dui
  21. #define        THIS_FFTU    dfft2du
  22. #define        THIS_SCAL    dscal2d
  23.  
  24. #endif
  25. #ifdef    SINGLE
  26.  
  27. typedef        float        this_half;
  28. typedef        float        this_type;
  29.  
  30. #define TOLERANCE   STOLERANCE
  31.  
  32. #define        THIS_SUM    ssum_
  33. #define        THIS_GEN    sgen_
  34. #define        THIS_FTU    sft2du_
  35.  
  36. #define        THIS_FFTUI    sfft2dui
  37. #define        THIS_FFTU    sfft2du
  38. #define        THIS_SCAL    sscal2d
  39.  
  40. #endif
  41.  
  42. /*                        */
  43.  
  44. this_half THIS_SUM();
  45.  
  46. void inimat_();
  47. void GetArguments();
  48. void get_values();
  49.  
  50. int size_x, ldx, size_y, n_trials;
  51. this_type *pa, *pb, *pRef, *pSave;
  52.  
  53. main(argc,argv)
  54. int argc;
  55. char *argv[];
  56. {
  57.     int i, j, k, n_total, is_wrong, nerr, inc;
  58.     double x, y, dx, dy, emax, sum, err, maxerr;
  59.     this_half    t;
  60.     this_type *ptr;
  61.  
  62. /* ******************************************************* */
  63.     GetArguments( argc, argv);
  64. /* ******************************************************* */
  65.  
  66.     srandom( (123*getpid()) | 0x01);
  67.  
  68.     for( ; n_trials > 0; n_trials --) { 
  69.     get_values();
  70.  
  71.     printf("\n");
  72.     printf( "%3d <%3d,%3d>", n_trials, size_x, size_y);
  73.     fflush(stdout);
  74.  
  75.     ldx = 2*((size_x + 2)/2);
  76.     n_total = (ldx*(size_y+1));
  77.     pa = (this_type *)malloc
  78.         ( (3 * n_total + size_x + 3 * size_y + 4 * FACTOR_SPACE) * sizeof(this_type));
  79.     if( !(pa) ) {
  80.         fprintf( stderr, "Could not allocate ... Exiting\n");
  81.         exit (-1);
  82.     }
  83.     pb = pa + n_total;
  84.     pRef = pb + n_total;
  85.     pSave = pRef + n_total;
  86.  
  87. /* *******************************************************
  88.     Initialize arrays
  89. ******************************************************* */
  90.     THIS_GEN(pRef, &n_total);
  91.     bcopy( pRef, pa, n_total*sizeof(this_type));
  92.     bcopy( pRef, pb, n_total*sizeof(this_type));
  93.  
  94. /* *******************************************************
  95.     NOT SO PACKED --- direct Fourier Transform
  96. ******************************************************* */
  97.         j = -1;
  98.     THIS_FTU( &j, &size_x, &size_y, pb, &ldx);
  99.     pSave = THIS_FFTUI( size_x, size_y, pSave);
  100.     bcopy( pRef, pa, n_total*sizeof(this_type));
  101.     is_wrong = THIS_FFTU( -1, size_x, size_y, pa, ldx, pSave);
  102.  
  103.     emax = TOLERANCE * TOLERANCE * (size_x * size_y);
  104.     nerr=0 ;
  105.     sum = 0.0;
  106.     err = 0.0;
  107.     for(i = 0; i < n_total ; i ++) {
  108.         sum += ( pa[i] * pa[i]);
  109.         x = pa[i] - pb[i];
  110.         t = x * x;
  111.               err += t;
  112.         if(t > emax) nerr++;
  113.     }
  114.     err = sqrt(err / (double)(size_x*size_y));
  115.     sum = sqrt(sum / (double)(size_x*size_y));
  116.     if( (err/sum) < TOLERANCE  ) {
  117.         printf(".");
  118.     }else {
  119.         printf(": avg:(%8.3g /%8.3g)= %8.3g >< ", err, sum, err/sum);
  120.         printf("\nDIRECT Fourier Transform : %d errors detected \n", nerr);
  121. /*
  122.         exit(-2);
  123. */
  124.     }
  125. /* *******************************************************
  126.     NOT SO PACKED --- inverse Fourier Transform
  127. ******************************************************* */
  128.     fflush( stdout);
  129.     is_wrong = THIS_FFTU( 1, size_x, size_y, pa, ldx, pSave);
  130.  
  131.     x = 0;
  132.     j = (size_y -1) / 2;
  133.     inc = 2 * ldx;
  134.  
  135.     t = 1. / (double)(size_x * size_y);
  136.     THIS_SCAL( size_x, size_y, t, pa, ldx);
  137.  
  138.     emax = TOLERANCE;
  139.     emax = emax * emax;
  140.  
  141.     sum = 0.;
  142.     err= 0.;
  143.     nerr= 0;
  144.     maxerr= 0.;
  145.     for ( j = 0; j < size_y ; j ++) { 
  146.       for(i = 0 ; i < size_x ; i ++) {
  147.         sum += (pRef[i+j*ldx] * pRef[i+j*ldx]) + (pRef[i+j*ldx] * pRef[i+j*ldx]);
  148.         x = pa[i+j*ldx] - pRef[i+j*ldx];
  149.         if( (t= x*x) > (emax))
  150.         nerr++;
  151.         err += t;
  152.         if( t > maxerr) maxerr = t;
  153.       }
  154.     }
  155.     err = sqrt(err / (double)(size_x*size_y));
  156.     sum = sqrt(sum / (double)(size_x*size_y));
  157.     maxerr = sqrt(maxerr);
  158.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g", 
  159.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  160.     if(nerr){
  161.         printf("\n%d errors detected \n", nerr);
  162.         exit(-2);
  163.     }
  164.     free(pa);
  165.     }
  166.     printf("\n");
  167.     return(0);
  168. }
  169.  
  170. int is_random;
  171.  
  172. void GetArguments( argc, argv)
  173. int argc;
  174. char *argv[];
  175. {
  176.     int i, j, k;
  177.     int nerror = 0;
  178.  
  179. #define ON    1
  180.  
  181.     n_trials = DEF_TRIALS;
  182.     is_random = 1;
  183.  
  184. /* ******************************************************* */
  185.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  186.     if( argv[i][0] == '-') {
  187.         switch ( argv[i][1]) {
  188.         case 'i' :
  189.         case 'I' :  
  190.                 is_random = 0;
  191.                 break;
  192.         default  :  nerror = ON;
  193.         }
  194.     }
  195.     else {
  196.         if( i+1 > argc)
  197.         nerror = ON;
  198.         else { 
  199.         n_trials = atoi( argv[i]);
  200.         }
  201.     }
  202.     }
  203.     if( nerror == ON) {
  204.     fprintf( stderr, 
  205.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  206.     exit(-1);
  207.     }
  208. }
  209.  
  210. void get_values()
  211. {
  212.   if( is_random){ 
  213.     size_x = random() % MAX_2D + 1;
  214.     size_y = random() % MAX_2D + 1;
  215.     if( !(n_trials % 5))
  216.     size_y = size_x;
  217.   } else{
  218.     printf( "Enter SIZE_X : ");
  219.     scanf("%d", &size_x);
  220.     printf( "Enter SIZE_Y : ");
  221.     scanf("%d", &size_y);
  222.   }
  223. }
  224.